Surface structure in simple liquid metals. An orbital free first principles study. 
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Molecular dynamics simulations of the liquid-vapour interfaces in simple sp-bonded liquid metals 
have been performed using first principles methods. Results are presented for liquid Li, Na, K, 
Rb, Cs, Mg, Ba, Al, Tl, and Si at thermodynamic conditions near their respective triple points, 
for samples of 2000 particles in a slab geometry. The longitudinal ionic density profiles exhibit 
a pronounced stratification extending several atomic diameters into the bulk, which is a feature 
aheady experimentally observed in liquid K, Ga, In, Sn and Ifg. The wavelength of the ionic oscil- 
lations shows a good scaling with the radii of the associated Wigner-Seitz spheres. The structural 
rearrangements at the interface are analyzed in terms of the transverse pair correlation function, 
the coordination number and the bond-angle distribution between nearest neighbors. The valence 
electronic density profile also shows (weaker) oscillations whose phase, with respect to those of the 
ionic profile, changes from opposite phase in the alkalis to almost in-phase for Si. 

PACS numbers; 61.25.Mv, 64.70.Fx, 71.15.Pd 
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I. INTRODUCTION 

The study of the structure of the free liquid sur- 
face has attracted much theoretical and experimental 
work^*^ with emphasis on the possible existence of liq- 
uid surface layers. Although it is by now well estab- 
lished that ionioi and dielectricii liquids exhibit just a 
smooth monotonic decay from the bulk liquid density to 
the bulk vapour density, liquid surface layering appears 
in liquid crystals^ and at the interface between a simple 
fluid and a hard waY^. Understandably, special atten- 
tion has been devoted to liquid metals since the early 
experiments suggested a liquid surface layered structure 
in Hgi; computer simulationSiSiifliiiiiSii^ and theoretical 
studies^'* predicted a significant structure, extending sev- 
eral atomic layers into the bulk liquid. These results have 
been recently corroborated by X-ray reflectivity experi- 
ments on liquid Hg, Ga, In, Sn, K and Nao.33 Ko.e?"^ 

It not yet settled whether surface layering is a fea- 
ture of all liquid metals or only characteristic of some. 
Rice and coworkersSiSiifi have performed Monte Carlo 
(MC) simulations based on density-dependent pair po- 
tentials obtained from pseudopotentials; their results 
suggest that surface layering is due to the coupling be- 
tween ionic and electronic density profiles (DP) and that 
the abrupt decay of the electron DP gives rise to an ef- 
fective wall against which the ions, behaving like hard- 
spheres, stack. Other workers^^ have suggested that the 
many body forces arising from the delocalized electrons, 
tend to increase the ionic surface density so that its coor- 
dination approaches that of the bulk. Recently, Chacon 
et afii proposed that surface layering may be a generic 
property of fluids at low temperature, so that the only 
requirement for an oscillatory DP is a low melting tem- 
perature relative to the critical temperature so as to avoid 
crystallization. 

This paper reports a study on the liquid-vapour in- 



terface of several simple sp-bonded liquid metals at ther- 
modynamic conditions near their respective triple points. 
The study has been carried out by using the orbital-free 
ab initio molecular dynamics (OF-AIMD) method, where 
the forces acting on the nuclei are computed from elec- 
tronic structure calculations, based on density functional 
theory (DFT)l^, which are performed as the MD trajec- 
tory is generated. Large samples and long MD simulation 
times are possible when interatomic pair potentials are 
used to describe the effective ion-ion interactions; how- 
ever at a liquid metal surface this approximation be- 
comes unreliable as the electron density sharply drops 
from its bulk value to zero outside the surface. In these 
circumstances it is vital to allow the forces on the atoms 
to respond to the electron density distribution in their 
vicinity. At present the best way of accomplishing this 
goal is by resorting to flrst-principles molecular dynamics 
techniques, where the electronic density, total energy and 
forces are obtained by using the Kohn-Sham (KS) formu- 
lation of the density functional theory (DFT).l8j However, 
the computational demands of these ab-initio methods, 
where KS orbitals are used to describe the electronic den- 
sity and to compute exactly the electronic kinetic en- 
ergy, grow very rapidly with system size, and their mem- 
ory requirement is also quite large. These considerations 
have posed important constrains on both the sizes of the 
systems studied so far, as well as the simulation times. 
Heretofore, only two KS-AIMD calculations have been 
performed on the liquid-vapour (LV) interfaces of liquid 
metals. Fabricius et alii have studied the LV interface 
in liquid silicon near melting by using 96 particles and a 
total simulation time of 30 ps. Likewise, Walker et 
have studied the LV surface in liquid sodium near melt- 
ing by a simulation which used 160 particles and lasted 
50 ps. In both studies, the size limitations materialized 
in small simulation slabs (~ 12.0 A thick for Si and 22.0 
A thick for Na) which may rise some reasonable ques- 
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tions about its capability to simulate a real macroscopic 
LV interface, as the small thickness makes plausible the 
existence of interactions between both sides of the slab. 

However, the aforementioned limitations can be partly 
overcome if the exact calculation of the electronic ki- 
netic energy is given up in favor of an approximate ki- 
netic energy functional of the electronic density. This is 
the philosophy behind the so-called orbital-free ab initio 
molecular dynamics (OF-AIMD) mcthodi^iSS, which by 
eliminating the orbitals of the KS formulation provides 
a simulation method where the number of variables de- 
scribing the electronic state is greatly reduced, enabling 
the study of larger samples and for longer simulation 
times. Although the OF-AIMD method uses an explicit 
approximation for the electronic kinetic energy functional 
of DPT, nonetheless it correctly treats the forces on the 
ions. We stress that, although less demanding than KS- 
AIMD, the OF-AIMD method is still much more costly 
than the use of pair potentials; however, recent calcu- 
lations have succeeded to use 2000 and 3000 particles 
to study the surface properties of several simple liquid 
metals and binary alloys, namely, Li, Na, Mg, Al, Si, 
Nao.aKo.T and Lio.4Nao.6. 

The OF-AIMD results for Na and Si, when compared 
with experimental data (as the surface tension) and with 
the in principle more accurate KS-AIMD results (as the 
density profiles) for the same systems, can be used as 
a validation test of the method, as far as the study of 
metallic surface properties is concerned. As we show be- 
low in detail, several magnitudes are very similar: the 
oscillating profiles, and in particular the wavelength of 
the oscillations, which is recovered exactly, the number 
of neighbors of a Si atom across the interface, and the 
surface tension of liquid Na, are all well reproduced by 
the OF-AIMD approach. Further confidence in the ca- 
pabilities of the OF-AIMD method in relation to surface 
properties is obtained from results for metallic clusters 
(finite systems where the surface is indeed utmost im- 
portant). For instance, a long standing, previously unex- 
plained, anomalous variation of the melting temperature 
of Na clusters with size^ has been for the first time repro- 
duced and rationalized in terms of surface geometry and 
stability^ using the OF-AIMD method with the same 
kinetic energy functional used in the present work. 

We do consider that reproduction of the surface prop- 
erties of finite systems is much more stringent a test than 
the reproduction of the surface properties of a semiinfi- 
nite solid surface. Nevertheless, in order to satisfy the 
requirements of one referee and to further reassess for 
the readers the capabilities of the OFAIMD method, we 
have performed preliminary studies of the properties of 
some open metallic solid surfaces, as the (110) surface 
of fee Al and the (1010) surface of hep Mg.^ In both 
cases, experimental measurements^iiS^ show that surface 
relaxation leads to a contraction of the first interlayer dis- 
tance, expansion of the second interlayer distance, con- 
traction of the third, and so on. Moreover, the thermal 
expansion coefficient in the case of Al (110) is negative. 



positive and positive for the first, second and third in- 
terdistances, while for Mg (1010) an oscillatory behavior 
is found starting with a negative thermal expansion co- 
efficient. KS-AIMD simulations for several temperatures 
were pcrformed^^ for Al (110) using 8 layers with 9 atoms 
in each one plus a vacuum of 8.5 A, reproducing the ex- 
perimental trends. Our OF-AIMD calculations, using 
the same simulation setup, also reproduce both the sign 
of the relaxations and their thermal behaviour. The- 
oretical calculations within the quasiharmonic approxi- 
mation, based on KS ab initio static calculations (not 
MD simulations) were also able to reproduce the ex- 
perimental behavior of Mg (1010). Our OF-AIMD data 
again reproduce qualitatively the experimental trends, 
with both oscillatory relaxations and thermal expansion 
coefficients. Further details and analysis of the data will 
be presented elsewhere. 

In our previous OF-AIMD studies of liquid metallic 
surfaces^ we first of all proved the feasibility of per- 
forming ab initio simulations for large systems, includ- 
ing one component metals and alloys, we showed that 
different ordering properties in the case of alloys lead to 
substantial differences in the partial and total density 
profiles, and studied the evolution of the relationship be- 
tween ionic and valence electron density profiles as the 
valence of the metal is increased in the series Li, Mg, 
Al, Si. Here we extend the number of systems studied, 
we analyze in detail the structure of the systems, both 
perpendicular and parallel to the interface, and complete 
the study of the electronic density profiles, stressing the 
evolution of these properties as the atomic valence is var- 
ied. 



II. THEORY 

A simple liquid metal is treated as a disordered array of 
N bare ions with valence Z, enclosed in a volume V, and 
interacting with N^, = NZ valence electrons through an 
electron- ion potential v{r). The total potential energy of 
the system can be written, within the Born-Oppenheimer 
approximation, as the sum of the direct ion-ion coulombic 
interaction energy and the ground state energy of the 
electronic system under the external potential created 
by the ions, Foxt(^ {Ri}) = Y^iLi , 

z<j - Rj\ 

where pg (f) is the ground state electronic density and Ri 
are the ionic positions. 

According to DFT, the ground state electronic den- 
sity, Pg{r), can be obtained by minimizing the energy 
fimctional E[p], which can be written 

E[p{f)] = T, [p] + Eh [p] + Sxc [p] + Soxt [p] (2) 
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where the terms represent, respectively, the electronic ki- 
netic energy, Tg [p\ , of a non- interacting system of density 
p(r), the classical electrostatic energy (Hartree term), 

the exchange-correlation energy, i?xc[p], for which we 
have adopted the local density approximation and finally 
the electron- ion interaction energy, i?oxt[p], where the 
electron-ion potential has been characterized by a local 
ionic pseudopotential which has been constructed within 
DFTia. 



E,Ap\ = j drp(f)V;xt(r), (4) 

In the KS-DFT method, Ts[p\ is calculated exactly 
by using single particle orbitals, which requires huge com- 
putational effort. This is alleviated in the OF-AIMD 
approacbi^iiSi^ by use of an explicit but approximate 
functional of the density for [p] . Proposed functionals 
consist of the von Weizsacker term, 

Tw[pm = IJ dr'|Vp(r)|Vp(r1, (5) 

plus further terms chosen in order to reproduce correctly 
some exactly known limits. Here, we have used an aver- 
age density modeUS, where Ts — Tw + Tp , 

Tp = yJ drp{r)''^-'0~m^ (6) 
fc(f) = {2k°pf j dsk{s)Wfj{2k°p\r-s\) 

k{r) = (Stt^)^/'^ Fermi wavevector for 

mean electron density pe = Ne/V ^ and wpix) is a weight 
function chosen so that both the linear response theory 
and Thomas-Fermi limits are correctly recovered. Fur- 
ther details are given in reference p9|. 

Another key ingredient of the energy functional is the 
the local ion pseudopotential, Wps(r), describing the ion- 
electron interaction. For each system, the Vps(r) has 
been constructed from first principles by fitting, within 
the same Ts[p] functional, the displaced valence elec- 
tronic density induced by an ion embedded in a metallic 
medium as obtained in a KS calculation which, more- 
over, also gives the corresponding core electronic density. 
Further details on the construction of the pseudopoten- 
tial are given in reference and we just note that the 
previous theoretical framework has already provided an 
accurate description of several static and dynamic prop- 
erties of bulk liquid Li, Mg, Al, Si, Na-Cs and Li-Na 



TABLE L Input data for the simple metals studied in this 
work, along with some simulation details. St is the ionic time 
step in ps, the cutoff energy, -Ecut, is given in Ryd, and A^conf 
is the total number of configurations. 



Metal 




T (K) Lo (A) 


a 


St 


Ecut 


iVconf 


Li 


0.0445 


470 


28.44 


1.95 


0.0060 


9.50 


18000 


Na 


0.0242 


373 


33.50 


2.20 


0.0025 


7.50 


24700 


K 


0.0127 


343 


44.99 


1.73 


0.0050 


5.25 


25100 


Rb 


0.0103 


315 


48.26 


1.73 


0.0065 


5.25 


22100 


Cs 


0.0083 


303 


51.77 


1.73 


0.0050 


4.74 


18300 


Mg 


0.0383 


953 


29.90 


1.95 


0.0010 


8.50 


22000 


Ba 


0.0146 


1003 


43.00 


1.73 


0.0040 


4.95 


21000 


Al 


0.0529 


943 


28.97 


1.56 


0.0010 


11.25 


20000 


Tl 


0.0332 


590 


32.53 


1.75 


0.0075 


10.50 


30000 


Si 


0.0555 


1740 


27.41 


1.75 


0.0035 


15.55 


20000 



III. RESULTS 

We have performed OF-AIMD simulations for the LV 
interfaces in the liquid metals Li, Na, K, Rb, Cs, Mg, Ba, 
Al, Tl and Si at thermodynamic conditions near their ex- 
perimental triple points. For each system we have con- 
sidered a slab consisting of 2000 ions in a supercell with 
two free surfaces normal to the z-axis. The dimensions 
of the slab were Lq ■ Lq ■ {Lz= a Lq), with Lq and a 
chosen so that the average number density of the slab co- 
incides with the experimental bulk ionic number density 
of the system at the same temperature; additional de- 
tails about the thermodynamic states are given in Table 
^ along with several simulation parameters. A further 8 
A of vacuum were added both above and below the slab. 
Therefore, we are dealing with liquid slabs which are wide 
enough to rule out interference effects between the two 
free surfaces and with supercells which are large enough 
to discard slab-slab interactions. Although the periodic 
boundary conditions require that a particle moving out 
of the cell in the z-direction reappears on the other side of 
the slab, we did not observe such event during the present 
simulations. Given the ionic positions at time t, the elec- 
tron density was expanded in plane waves and the energy 
functional was minimized with respect to the plane wave 
coefficients yielding the ground state electronic density, 
energy, and the forces on the ions, and therefrom the 
ionic positions and velocities were updated according to 
Newton equations, i.e., the simulations are performed in 
the NVE ensemble. For all systems equilibration runs 
were previously performed for a range between 2000-4000 
configurations, depending on the system. Therefrom, the 
Aconf ensuing configurations were those used in the eval- 
uation of the slab's physical properties. 

During the simulations each slab contracted or ex- 
panded and the average ionic density varied in response 
to the condition of zero external pressure so that the 
ionic density in the central region of slab changed by an 
amount ranging from -1.5% in Tl to « 20% in K and 
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FIG. 1: Electronic (dotted line) and ionic (full line) density 
profiles normal to the liquid-vapour interface in the Li, Na an 
K liquid metals. The densities are plotted relative to their re- 
spective bulk values. The dashed and dot-dashed lines are the 
a;-transverse (displaced by -0.25) and j/-transverse (displaced 
by -0.5) ionic density profiles. 



Si. 



A. Ionic density profiles. 

The longitudinal ionic DP were computed from a his- 
togram of particle positions relative to the slab's center 
of mass, with the profiles from both halves of the slab 
being averaged, and the obtained results are shown in 
figures ^121 All systems show a stratification for at least 
four layers into the bulk liquid, a structural feature that 
has already been observed experimentally for the LV in- 
terface of Hg, Ga, In, Sn, K and NasaKgr^^. 

The wavelength. A, of the ionic oscillations is given in 
Table Inl and its values scale linearly with the radii of the 
Wigner-Seitz spheres, Rws , while it shows no definite re- 
lationship with electronic parameters, like the radii per 
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FIG. 2: Same as the previous figure but for the Rb and Cs 
liquid- vapour interfaces. 



electron, (see figure El . This result suggests that the 
ionic oscillations are not a consequence of Friedel oscilla- 
tions in the electronic density, but on the contrary they 
are due to atomic stacking against the surface. Paren- 
thetically, we mention that the present OF-AIMD values 
of A for Na (3.0 A) and Si (2.5 A) exactly coincide with 
the values derived in KS-AIMD type calculationspii^i^ 

Another noticeable feature of the ionic DPs refers to 
the relative amplitudes of the first and second oscillations 
near the LV interface. Within the group of the alkali met- 
als, we observe that the amplitude of the first oscillation 
is smaller than that of the second oscillation; this feature 
is more marked in Li and is much smaller for the other 
alkalis. Conversely, for Mg we obtain a first oscillation 
whose amplitude is slightly larger than that of the sec- 
ond albeit this effect is reversed in Ba. However, for Al, 
Tl and Si we obtain a first oscillation with an amplitude 
clearly larger than that of the second oscillation. Indeed, 
the present results suggest that as the valence of the sys- 
tem increases so does the ratio between the amplitudes of 
the first and second oscillations in the longitudinal ionic 
DP. 

The OF-AIMD results for the alkalis closely agreee 
with previous MC results of Rice et a&2ii£ for Na, K, 
Rb and Cs, where the first oscillation had a shorter am- 
plitude than the second. Furthermore, their MC results 
for Al, Ga, In and TiM showed a first oscillation with 
a larger amplitude than the second, which again coin- 
cides with the OF-AIMD results for Al and Tl, and most 
importantly, this trend is also visible in the reported ex- 
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FIG. 3: Same as the previous figure but for tfie Mg and Ba 
liquid- vapour interfaces. Notice tfie different a;-ajds scales. 
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FIG. 4: Same as the previous figure but for the Al and Tl 
liquid- vapour interfaces. Notice the different a;-axis scales. 
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FIG. 5: Same as the previous figure but for the Si liquid- 
vapour interfaces. 



TABLE II: Information about several properties of the slabs. 
p is the ionic number density in a wide central region of the 
slab, Ap is its variation with respect to the input ionic density, 
A is the wavelength of the longitudinal ionic oscillations, and 
a denotes the pseudoatom size (see text). 

Metal p (A-^) Ap (%) A (A) a (A) a/A 



Li 


0.0478 


7.6 


2.5 


1.601 


0.640 


Na 


0.0272 


12.0 


3.0 


1.859 


0.620 


K 


0.0152 


19.4 


3.7 


2.309 


0.624 


Rb 


0.0118 


14.4 


3.9 


2.458 


0.630 


Cs 


0.0089 


7.2 


4.3 


2.710 


0.630 


Mg 


0.0404 


5.5 


2.6 


1.500 


0.577 


Ba 


0.0162 


11.1 


3.6 


2.120 


0.589 


Al 


0.0570 


7.2 


2.35 


1.288 


0.548 


Tl 


0.0328 


-1.5 


2.90 


1.274 


0.440 


Si 


0.0666 


20.0 


2.5 


1.132 


0.453 



perimental longitudinal ionic DP for Hg, Ga and Ini^. 
Likewise, we note that both the present OF-AIMD and 
the MC results, have yielded a first oscillation which is 
more marked in Tl than Al. On the other hand, the re- 
cent KS-AIMD simulatioriii. for the LV interface in Na, 
has also produced an oscillatory longitudinal ionic DP 
with a first oscillation having an amplitude slightly larger 
than that of the second. 

Besides the correlation found between the amplitudes 
of the first and second oscillations and the valence of the 
system, we have also unveiled another correlation con- 
cerning the decaying tail of the ionic DPs in the LV in- 
terface region Specifically, when the ionic DPs are scaled 
in terms of their respective wavelength. A, is it found that 
the decaying tail of p{z/X) is virtually the same for all 
the systems with the same valence, Z, and moreover this 
tail becomes steeper with increasing Z . 

We have also checked that the stratification of the ionic 
DPs is not an artifact induced by the finite size of the 
simulation box. Therefore, we have computed the trans- 
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TABLE III: Details of the layers used for computing the trans- 
verse pair correlation functions, ^zql and /S.zsl refer to the 
positions, with respect to the center of mass, of the outermost 
and second layers, respectively. Apoi and ApsL are the per- 
cent variations, with respect to the slab's bulk value, of the 
ionic number densities in the outermost and second layers, 
respectively. 



Metal AzoL (A) Azsl (A) ApoL 


ApSL 


Li 


23.6-25.7 


21.1-23.6 


-11.3 


-2.3 


Na 


29.9-32.4 


26.9-29.9 


-6.9 


0.0 


K 


29.3-32.4 


25.6-29.3 


-8.1 


-0.5 


Rb 


32.8-36.1 


28.9-32.8 


-7.2 


0.0 


Cs 


37.9-41.7 


33.6-37.9 


-9.4 


0.0 


Mg 


25.2-27.6 


22.6-25.2 


-6.5 


-1.0 


Ba 


30.2-33.4 


26.6-30.2 


-7.7 


-0.6 


Al 


18.7-20.8 


16.35-18.7 


+4.0 


-1.0 


Tl 


26.0-28.5 


23.1-26.0 


-1-9.0 


-1.2 


Si 


17.5-19.8 


15.0-17.5 


+0.6 


-1.5 




FIG. 6: Wavelength of the longitudinal ionic oscillations as 
a function of (a) the radius of the respective Wigner-Seitz 
spheres, Rws (b) the radius of a sphere which on average 
contains one electron, r^. 



verse ionic DPs which, as shown in figures [Tl5l are more 
or less uniform, albeit with some noise which is always 
substantially smaller than the amplitudes of the oscilla- 
tions in the corresponding longitudinal DP. As a further 
check on the reliability of the present calculations, we 
have also calculated the corresponding bulk pair distri- 
bution functions, g(r), which have been evaluated within 
a 30.0 A wide central section of the corresponding slab. 
The obtained results are depicted in figures [7[ along 
with their experimental counterparts^^. The small mis- 
match observed in Mg, Ba and Al may be ascribed to 
the aformentioned increase in the average ionic density 
in the central part of the slab. 



B. Atomic structure of the layers. 




FIG. 7: (Color online) Transverse pair correlation functions 
for selected layers of the Li and Na slabs, namely from the 
bulk (thin blue line) and from the outermost layer (thick line). 
The full circles stand for the bulk experimental data. 




The planar LV interfaces considered in this work are in- 
homogeneous systems along the z direction, but isotropic 
and homogeneous on the xy planes; consequently the 



FIG. 8: (Color online) Same as the previous figure but for 
liquid Mg and Ba. 
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FIG. 9: (Color online) Same as the previous figure but for 
liquid Al, Tl and Si. In Tl the dashed red line is the grir) of 
the second layer. 



two-body distribution functions, which provide informa- 
tion about the atomic structure, depend on the z coor- 
dinates of the two particles and on the transverse dis- 
tance between them, i.e. 17(^1,^2) — g{zi, Z2, R12) where 
i?i2 = {x2 — xiY + (2/2 — yiY- The three body distri- 
bution functions, as for instance the bond- angle distribu- 
tion, also display similar symmetry properties. 

However, such a detailed description of the structure 
is difficult to obtain from the simulations, so below we 
consider several averaged magnitudes in which the aver- 
age is taken over a given layer or slice of the simulation 
slab. The first (or outermost) layer comprises the re- 
gion from the outermost minimum of the ionic DP to 
its inflection point in the decaying tail. The other lay- 
ers are located between consecutive minima of the ionic 
oscillations. For a given metal, all the layers have the 
same width (excepting the outermost one which is al- 
ways slightly narrower) and its specific values are given 
in Table ITTll Furthermore, in order to achieve meaning- 
ful comparisons we have also considered a layer with the 



same width and located in the central bulk region, which 
will be termed the "center layer" . 

A first magnitude to analize is the average ionic den- 
sity within each layer. Another one is the transverse pair 
correlation function 5t(?'), which accounts for the prob- 
ability of finding two particles separated by a distance r 
when both particles are inside the layer. Then, we will 
consider the number of nearest neighbors (NN) and the 
distribution of angles between triplets formed by one par- 
ticle and its two NN. Notice that in order to describe the 
in-plane structure of a layer the two previous magnitudes 
should be evaluated for those neighbours standing inside 
the layer. Moroever, the particles in a layer have addi- 
tional neighbors in adjacent layers and therefore a more 
thorough description of the structure will require to ac- 
count for them in evaluating the total number of NN and 
the distribution of angles. 

1. Average ionic density. 

The calculation of the mean ionic number density 
within each slice has revealed that deviations from the 
bulk value are significant only in the last two slices. At 
the outermost slice the relative changes in the ionic num- 
ber density range from w -11 % (Li) to 9 % (Tl) whereas 
in the previous slice the corresponding values are sub- 
stantially smaller (see Table 1111)1 . Note that the average 
ionic density in the outermost layer decreases in all the 
alkalis and alkaline-earths, and increases for Tl, Al and 
Si. 



2. Transverse pair correlation functions. 

In all the systems studied, excepting Tl, the change in 
grir) from the bulk to the surface occurs rather abruptly 
at the outermost layer, because at the second layer its 
griT) already coincides with that of the center layer, 
which in turn is practically identical to the bulk g{r). 
As for Tl, its LV interface is the most structured one, 
with large oscillations which may indicate an important 
influence of any given layer on the properties of its sur- 
rounding layers. In fact, the change in grif) from the 
bulk to the surface is gradual in the case of Tl. 

The changes undergone by the average ionic density 
(Table lllljl are mirrored by deviations of the associated 
grir) with respect to the bulk g{r), as evinced by figures 
I7I9I which show the transverse pair correlation function 
grir) at the outer layer. For both alkalis and alkali- 
earths, the outermost gxif) is displaced towards greater 
r-values, in clear correlation with a decreased average 
ionic density. Conversely, the outermost griT) for Al, Tl 
and Si practically preserve the same main peak's position 
as in the bulk while slightly increasing its height. Notice 
that in these latter three systems, the ionic number den- 
sity at the outermost layer did increase with respect to 
the slab's bulk value. 
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TABLE IV: Number of neighbors within a layer n2d for the 
systems studied. 
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0.1 
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Pi 1 

0.1 


Cs 


4.7 


5.2 


5.2 


5.2 


Mg 


4.9 


5.1 


5.1 


5.1 


Ba 


4.8 


5.3 


5.3 


5.3 


Al 


5.1 


4.9 


4.8 


4.8 


Tl 


5.5 


5.5 


5.2 


5.1 


Si 


4.1 


3.9 


4.1 


4.0 



Conceivably, it might appear that for the alkalis and 
alkaline-earths those changes in grif) just amount to an 
expansion of the system, leading to an increased NN's 
distance. In the case of Al and Tl the maximum of grif) 
does not change position, but the first minimum moves 
somewhat to the right. The most interesting situation 
occurs for Si, where apparently some atoms move from 
the position of a small bump just after the main peak 
of g{r) to a distance slightly smaller than the position 
of the second maximum, yielding a double-second-peak 
structure in the grir) of the outermost layer. All these 
changes induced by the layering of the interface will be 
thoroughly analyzed in the following section. 



3. In-plane neighbors. 

The number of NN, or coordination number (CN), is 
usually defined as the average number of neighbors within 
a distance Tj^axi identified as the position of either the 
first minimum of the pair correlation function 5t(^), or 
that of the radial distribution function Grir) (which for 
these quasi-bidimensional layers is proportional to rgx (?') 
and the average ionic density) of the layer considered. 
This latter criterion has been used to calculate the re- 
sults of Table HVl although similar trends are found when 
using the former one. In the alkalis and alkaline-earths 
there is a competition between the increase of r^^^^ at 
the outermost layer and the corresponding decrease of 
the average ionic density. In all cases the latter factor 
is stronger and we find that the CN, n2d changes from 
about 5.2 in the center layer to about 4.7 in the outer- 
most one. In Al and Tl both the average ionic density 
and rf-^i^^ increase leading to a somewhat increased CN 
from the center to the surface. In Si both factors change 
very little and we find 4 NN (in-plane) both for the center 
layer and for the outermost one. 

Notwithstanding the validity of these results it must be 
recognized that there is a certain degree of arbitrariness 
in the definition of when a particle is a "nearest" neighbor 




FIG. 10: (Color online) Decomposition of the pair distribu- 
tion functions for Li in terms of those corresponding to the 
first, second, • ■ • ,up to the 20th neighbor. Thick lines: out- 
ermost layer, thin blue lines: center layer. Also shown are the 
associated pair correlation functions, grir) (scaled so as to fit 
into the graph). 



of another. Some years ago McGreevy and coworkersSS 
outlined a proposal for a less arbitrary determination of 
the CN in computer simulation studies. Their idea is to 
rewrite the radial distribution function as the sum of the 
partial functions, Gt{t) = X^i^i^rC*^)! where G\^{r) is 
the pair distribution function corresponding to the i-th 
neighbor, which indicates the distribution of distances 
from one particle to its i-th neighbor. By analyzing the 
way the different neighbors are distributed, it is possible 
to discern whether they are inside the first peak, outside 
it or in between, and therefore it provides a more precise 
picture of the meaning of the CN. In the present context 
this idea is even more useful, since by looking at the 
change of G'ipir) across the LV interface we can better 
describe the rearrangement of the particles induced by 
the layering of the interface. 

In figures [T0I12I we show this decomposition of Gt{t) 
for Li (which stands as a representative example of the 
alkalis and alkaline-earths, because all of them show the 
same behavior), Al and Si. Note that in all cases the 
amplitudes of G\<{r) follow closely the shape of the pair 
correlation function, grif), of the layer considered. 

Let us consider first the case of Li at the center layer. 
It is found that the neighbors up to the fourth are well be- 
low the first peak oi gxir), and the fifth neighbor, despite 
being more spread out, can reasonably be assigned to the 
set of "nearest" neighbors. The sixth one, however, is 
clearly distributed between the first and second peaks of 
grif)- The second peak is spawnned by part of the sixth 
neighbor, the neighbors number 7 to 15 and again partly 
the 16th one. From this analysis one could assign 5.4 and 
10.1 neighbors to the first and second coordination shells 
respectively. In the case of the outermost layer we find 4 
neighbors well inside the first peak of griT), whereas the 
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r(A) 

FIG. 11: (Color online) Same as the previous figure but for 
Al. 




r(A) 

FIG. 12: (Color online) Same as the previous figure but for 
Si. 

fifth and sixth neighbors are very spread out, both being 
shared by the first and second peaks of grir). The rest of 
the second peak is due to neighbors number 7 to 15 and 
part of the 16th one. We can therefore ascribe values of 
5 and 10.5 to the first and second CNs. In order to sep- 
arate out the effects of the decrease of the average ionic 
density, the comparison between the G'^(r) of the center 
layer and that of the outermost one is made for Li by 
rescaUng the distances by a, which denotes the position 
of the maximum oi gxir). In this way we observe that the 
expansion is not uniform: the peak positions of (r) to 
Gip (r) remain unchanged (in units of a) but from the 4th 
neighbor onwards there is more than uniform expansion 
and already the 6th neighbor in the outermost layer is 
located at the position of the 7th one in the bulk, while 
the 12th neighbor in the outermost layer is located at the 
position of the 14th one in the center. 

For liquid Al the rearrangement of atoms at the inter- 



face is very small. The first five neighbors hardly change, 
whereas starting at the sixth neighbor a small expansion 
occurs, followed again by contraction so that the posi- 
tions of the 15th neighbors again coincide. Analyzing the 
distribution of atoms under the first and second peaks, 
the CNs for the first and second shell would be 5.0 and 
9.5 in the center and 5.2 and 10.3 at the surface. 

In Si we find that the small bump after the main 
peak of the bulk g(r) is due to the particular distribu- 
tion of the fifth and sixth neighbors, the first four be- 
ing clearly located under the main peak, which leads to 
a CN of 4.0 atoms, both in the center and at the sur- 
face. When one comes to the outermost layer the heights 
of G^ (r) and (r) are depleted, whereas the following 
neighbors up to the 16th attain very similar distribu- 
tions although of course at increasing distances, yielding 
this way the double-peak structure of the second peak 
of gxir)- Therefore the previous interpretation of some 
atoms moving from the bump to the second peak is not 
correct; the position of the fifth neighbor is exactly the 
same in the surface and in the bulk, and starting from 
the sixth neighbor a weak expansion does occur, but it 
amounts to one atom only at large distances: the 18th 
surface neighbor is located at the position of the 19th 
bulk neighbor. 



4- In-plane bond-angle distribution functions. 

The usual way to calculate the bond-angle distribution 
function is to choose an atom and compute the angle be- 
tween the lines joining this atom and two of its NN, which 
would mean that they are separated less than r^^^ from 
the first one. However, following the ideas of the previous 
section we compute the angle between the lines joining 
the central atom and any two atoms taken from its first i 
neighbors. It is interesting to analyze the evolution of the 
bond-angle distributions, Pi{9), as more NN are included 
by starting at i = 2 and up to the integer closest to the 
CN. Figure E| shows the obtained results for liquid Li 
and Si. 

Starting with Li, we observe that in the center layer 
the distribution of cosines of angles peaks near 0.5, —0.5 
and —1 (angles 60, 120 and 180 degrees) suggesting that 
the local environment is mainly hexagonal, even though 
the CN is around 5 (in a penthagonal environment the 
peaks would occur at 0.309 and —0.809 very far from the 
observed ones). As we increase the number of neighbors 
included in the calculation the peak around 60 degrees 
moves towards smaller angles, whereas the one near 120 
degrees becomes fuzzier while hardly changing its posi- 
tion. The picture in the outermost layer is rather similar 
though the peak near 120 degrees is reduced to a shoulder 
in favor of an enhanced peak at 180 degrees. 

The case of Si, with a more open structure, shows a 
much stronger variation in the bond-angle distributions 
as more neighbors are included in the calculation. Start- 
ing with the center layer, we observe a clear peak near 
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TABLE V: Values of the z-dependent coordination number, 
n{z) at diferent ^-values along the slab. n{zB) is the average 
value in a wide region around the center of the slab, zqm 
and zsM are the positions (in A) of the outermost and second 
maximum respectively, z^i is the distance (in A) between the 
outermost maximum and the inflection point at the decaying 
ionic density profile, and n{zd) is defined in the text. 

Metal n{zB) zqm n{zoM) zsm n{zsM) Zd n{zd) 



Li 


11.85 


24.7 


7.9 


22.4 


11.3 


1.05 


7.8 


Na 


11.60 


31.1 


7.9 


28.3 


11.2 


1.30 


7.6 


K 


11.90 


30.9 


7.9 


27.4 


11.5 


1.45 


7.8 


Rb 


12.00 


34.5 


8.1 


30.9 


11.7 


1.55 


7.8 


Cs 


12.00 


39.9 


8.0 


35.8 


11.5 


1.80 


7.8 


Mg 


11.50 


26.4 


7.9 


23.9 


11.3 


1.25 


7.7 


Ba 


12.00 


31.8 


8.1 


28.4 


11.6 


1.60 


7.9 


Al 


11.30 


18.7 


7.9 


16.35 


11.2 


1.05 


7.6 


Tl 


11.90 


27.5 


8.2 


24.6 


11.9 


1.0 


7.6 


Si 


5.9 


18.8 


4.5 


16.3 


5.8 


1.15 


4.3 



FIG. 13: (Color online) Bond-angle distribution functions for 
Li and Si at different layers of the interface. The number of 
neighbors included in the calculation increases from 2 towards 
the coordination number (5 for Li, 4 for Si) as the position of 
the peak near cos 6 = 0.5 moves to the right. The vertical lines 
indicate the cosines of the tetrahedral angle (109.5 degrees) 
and of 90 degrees. 



60 degrees, which moves towards smaller angles if more 
neighbors are included. However, if only the first two 
neighbors are included a faint feature also appears near 
the tetrahedral angle, when including up to the third 
neighbor, this features moves to an angle near 90 degrees, 
and if we consider all the neighbors up to the fourth, the 
feature completely disappears. The variations brought 
about at the surface layer are similar as those observed 
in Li: we find an enhanced peak at 180 degrees, while 
the small features observed at the center layer for inter- 
mediate angles show up as shoulders, which also move 
their position from thetrahedral to 90 degrees and to dis- 
appearance when including 2, 3 or 4 neighbors in the 
computation. 



5. Total number of neighbors. 

To analyze the local 3-dimensional structure of the 
atoms at the different layers, we have computed the corre- 
sponding z-dependent CN, n(z), by counting the average 
number of neighbors (whatever be their position) within 
a distance r^^^ taken as the position of the first mini- 
mum of the bulk pair distribution function (proportional 
to r^g{r)). Our results, which are summarized in Table 
Ivl show that for most of the slab n{z) remains invari- 
able and close to that of the bulk system and only very 



close to the interface, namely starting around the second 
outer maximum, n{z) begins to decrease. Aside from liq- 
uid Si, all the other systems have an average bulk CN, 
n{zB) ~ 12, which is reduced by « 1/3 at their respective 
outer maximum, namely n{zoM) ~ 8. Liquid Si has a 
noticeably smaller n{zB) ~ 6, decreasing to « 4.5 at the 
surface, because of some remmants of covalent bonding 
in the liquid, which induce an open structure with an ex- 
perimental value of around S.'S^i- 6.4^ neighbours. For 
comparison, we note that the KS-AIMD calculations of 
Fabricius et aiii for liquid Si yielded n{zB) ~ 6.4 which 
reduced to 5.3 at the outermost maximum. Indeed, this 
decrease in CN is a well known feature at surfaces, pro- 
duced by the lack of atoms outside the interface. In the 
next section we will discuss if (and how) the remaining 
neighbors redistribute in this outer part of the system. 

As for the structural rearrangements induced by the in- 
terface, Fabricius et have proposed to quantify those 
changes by comparison with an ideally terminated sur- 
face obtained by cutting abruptly the slab in the central 
region, i.e. at z = 0. Then, we evaluate n{z) at a dis- 
tance Zd which is approximately the distance between 
the outermost maximum and the inflection point in the 
respective ionic DP. The obtained values are shown in 
Table Ivl and they are slightly smaller than those at the 
outer maximum which suggests that the surface struc- 
tural rearrangement induces some minor increase of the 
CN of an ideally terminated surface. 

6. Total bond-angle distribution functions. 

Again, instead of the usual way of defining a NN as one 
separated by at most r'^^^ from the central one, we com- 
pute the bond angles between a particle and two of its 
first i neighbors, with emphasis on i being either 2 or the 
integer value closest to the (3-dimensional) CN. Figure 
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depicts the results obtained for the bond-angle distri- 
bution functions in Li and Si. Specifically, we show the 
results obtained by counting only the first two neighbors 
in the calculation as well as those obtained by consider- 
ing the corresponding CN, which amounts to 8 neighbors 
for Li at the outermost layer and 12 at the inner ones, 
and 4 neighbors for Si at the surface, and 6 at the second 
and center layers. 

In both systems we find that when the calculation in- 
cludes all the NN (the full CN), then the bond- angle dis- 
tribution at the surface closely resembles that at the inner 
layers, including the bulk; however, we stress that this 
result relies on taking the CN associated to each layer, 
otherwise the distributions would peak at somewhat dif- 
ferent angles. The more probable angles are somewhat 
smaller than those reflecting a perfect icosahedral envi- 
ronment in Li, and around 60 degrees and a somewhat 
smaller angle than the tetrahedral one in Si. 

On the other hand, if the calculation includes only the 
first two neighbors, then the obtained bond-angle distri- 
bution describes a situation where the disorder induced 
by the distance is somehow at its minimum level, thus re- 
flecting the "most local" bond distribution. The results 
shown on the left panels of figure El provide evidence 
that in the case of Si, the surface has practically no in- 
fluence on the bond-angle distribution which now shows 
a peak around 60 degrees and another one at exactly the 
tetrahedral angle. For liquid Li we find peaks close to the 
icosahedral angles in the inner layers, but at the surface 
we observe a redistribution of bonds with an increased 
number of bonds at around 60 degrees along with an im- 
portant depletion at larger angles. However, the peak's 
positions do not change appreciably and remain around 
the icosahedral angles. 



C. Electronic structure. 

Figures FlISI show the calculated self-consistent valence 
electronic DP which also exhibit some oscillations al- 
though with much smaller amplitudes than those of the 
associated longitudinal ionic DP. However, the most re- 
markable feature refers to the relative phases between the 
longitudinal ionic and valence electronic DP. According 
to Figures FlISI all situations are allowed with the relative 
phase evolving from an opposite one, as it happens for 
all the alkalis, to being almost in-phase for Si. 

This is an unexpected feature because the scant pre- 
vious studies on the longitudinal ionic and the valence 
electronic DP^''^^ had always yielded an opposite phase. 
Specifically, the Monte Carlo calculations of Rice and 
coworkers^ for the LV interface in liquid Na, Cs and Ga 
produced longitudinal ionic and valence electronic DP 
standing in nearly perfect opposite phase. This feature 
was rationalized in terms of a competition between the 
valence electronic kinetic energy contribution, which gets 
smaller values by damping the oscillations, and the in- 
teraction term between electrons and ions, (mainly its 



i=2 i= coordination 



0.015 
0.01 
0.005 


0.006 
0.003 



- 




1 


,11,11,11,1 




I tvi\ - 

1 1 1 1 1 1 pTY 


i 






i 






\ 





1 






% 60 90 120150 30 60 90 120150180 



FIG. 14: (Color online) Total bond-angle distribution func- 
tions for Li and Si. Thick line: outermost layer. Dashed red 
line: second layer. Thin blue line: center layer. The number 
of neighbors included in the calculation is 2 for the left panels 
and the coordination number for the right panels. The arrows 
indicate the icosahedral angles (63.4 and 116.6 degrees) for Li 
and the tetrahedral angle (109.4 degrees) for Si. 



coulombic part), which being attractive takes smaller val- 
ues for in-phase oscillations. Furthermore, the OF-AIMD 
simulations by Jesson and Madden'^'^ for the structure of 
the liquid-solid interface in Al also gave a longitudinal 
valence electronic DP standing in a opposite phase with 
the ionic one in both the crystalline and liquid phases. 
This behaviour was now explained in terms of the interac- 
tion, embodied in the pseudopotential, between valence 
and core electrons, which tends to exclude the valence 
electrons from the ionic positions. 

In fact, the assumption of an opposite phase between 
the longitudinal ionic and valence electronic DP has been 
widely assumed; however a closer scrutiny, based on the 
present calculations, reveals some weak points on the pre- 
vious explanations. First, our results for the valence elec- 
tronic kinetic energy show that it takes very small values 
(less that 5% for all systems in this work) in comparison 
with those of the valence electron-ion interaction term; 
therefore the idea of assigning to the valence electronic 
kinetic energy any relevant role in establishing the phase 
of the valence electronic DP 'oscillations appears base- 
less. Second, all the systems considered in this work, 
have s-type valence electrons which in some cases (Al, 
Tl and Si) are also joined by p-type electrons. Taking 
into account that the s-type valence electronic density 
attains a maximum value at the ionic positions, it now 
appears rather surprising to find systems, specially those 
with s-type valence electrons only, in which both profiles 
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FIG. 15: Valence pseudoatomic densities for K (dotted line), 
Mg (dashed line) and Si (full line); the inset depicts its shape 
(multiplied by a factor of ten) for larger distances. The open 
circles represent the analytical model proposed in the text. 



were not in-phase. Consequently, we have undertaken 
several tests in order to ascertain the mechanisms under- 
lying the phase-shift between the longitudinal ionic and 
valence electronic DP. 

First, we have constructed a valence electronic DP in 
the slab by performing a linear superposition of the va- 
lence pseudoatomic density calculated within the process 
of constructing the local ionic pseudopotential. We re- 
call that the valence pseudoatomic density, which was 
obtained by means of a KS-type calculation, is a spheri- 
cally symmetric function and in order to analize its role 
in setting the shape of the valence electronic DP, we have 
first integrated over its x and y-dependencies; the result- 
ing function 7ie(z), which is the one effectively used in 
the calculation of the valence electronic DP in the slab, 
is plotted for some representative systems, in figure 1151 
For all systems, the corresponding ne{z) shows a decay- 
ing behavior whose width along with the associated weak 
Friedel oscillations may be identified as its basic charac- 
teristic features. The calculation of the valence electronic 
DP in the slab has proceeded as follows: by taking the 
OF-AIMD generated ionic positions we have located at 
each ionic position the previous ne{z) and we have carried 
out a configurational average. This approach is consis- 
tent with a linear response treatment of the valence elec- 
tron density and therefore lacks any competition between 
kinetic and coulombic effects. The resulting valence elec- 
tronic density profiles are compared with the OF-AIMD 
self-consistent ones for K, Mg and Si in figure^] which 
shows that both profiles are remarkably similar to each 
other. In particular, the phase of the valence electronic 
density oscillations is preserved, which suggests that the 
phase-shift must be determined by some feature embod- 
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FIG. 16: Valence electronic density profiles (relative to the 
bulk average value) for the K, Mg and Si liquid-vapour inter- 
faces. The continuous line is the OF-AIMD result, the dashed 
line represents a linear superposition of displaced densities 
and the grey circles are a linear superposition of the model 
densities, without Friedel oscillations. 



ied in the valence pseudoatomic density. 

In order to analize the different roles played by both 
features in modulating the phase-shift between the lon- 
gitudinal ionic and valence electronic DP, we have fitted 
the ne(z) to a model that shows no Friedel oscillations 
but otherwise accurately reproduces its shape. We found 
that a rather good fit can be obtained for a model den- 
sity of the normalized form exp[— jz/crp], which includes 
just one parameter, a. The specific values for the all 
the different systems are given in Table ^ which shows 
a range of variation from a = 2.71 A in Cs to cr = 1.132 
A for Si. Furthermore, a comparison between the ne{z) 
and the above fitting model, is depicted in Figure [T31 for 
K, Mg and Si. Besides the absence of the Friedel oscil- 
lations, it is observed that the model provides a rather 
good description of n^i^z) for all the systems in this work. 
Now, by using the fitted model densities we have again 



performed its superposition according to the OF-AIMD 
generated ionic positions and the obtained valence elec- 
tron profiles are shown in figure 1161 for K, Mg and Si. 
Once again we observe that the corresponding valence 
electronic DP is virtually indistinguishable, for all sys- 
tems, from that obtained by the superposition of valence 
pseudoatomic densities. More interestingly, the phase of 
the oscillations is preserved, and this fact suggests that 
the Friedel oscillations of the valence pseudodensity are 
irrelevant for this question. Therefore, the reason for the 
different phase-shifts between ionic and valence electron 
oscillations must be found in the width of the valence 
pseudoatomic density (quantified by the parameter a) as 
compared to the distance between layers in the profile 
(quantified by the wavelength of the longitudinal ionic 
oscillations A). The ratio cr/A takes the values 0.62, 0.58 
and 0.45 for K, Mg and Si respectively, which moreover 
correlates with a decreasing phase difference between 
their associated longitudinal ionic and valence electronic 
DPs. Table ITTl lists the obtained values of a/X for all the 
systems considered in this work and three groups can be 
discerned: (i) the alkalis where 0.62 < u/A < 0.64, (ii) 
Mg, Ba and Al with 0.55 < a/X < 0.59, and (iii) Tl and 
Si where 0.44 < a/X < 0.47. Within each group, there is 
a similar phase difference between the corresponding lon- 
gitudinal ionic and valence electronic DPs as epitomized 
by K, Mg and Si which may be considered as represen- 
tatives for each group. 

As a further check on the previous argument, we have 
performed another test for the K, Mg and Si slabs. By 
taking their respective OF-AIMD generated ionic posi- 
tions, we have again carried out a superposition of model 
densities with different widths, namely tr/A « 0.45, 0.575 
and 0.625, which may be considered as representative val- 
ues of the previous three groups. The calculated valence 
electronic DPs are depicted in figureElwhich shows that, 
in the three systems, the resulting valence electronic DPs 
move from being nearly in-phase (when a/X ~ 0.45) to- 
wards nearly opposite phase (when a/X « 0.625) with re- 
spect to the corresponding longitudinal ionic DP. These 
results clearly showcase the crucial role played by the 
ratio a/X in establishing the phase-shifts between the 
longitudinal ionic and valence electronic DPs. Moreover, 
the previous results provide a rationale for the aforemen- 
tioned phase-shifts between the ionic and valence elec- 
tronic DPs obtained for liquid Al in the present cal- 
culations and the OF-AIMD simulations of Jesson and 
Madden^. Specifically, our OF-AIMD results for the 
liquid-vapor interface of Al yielded A — 2.35 A, (with a 
ratio a/XK, 0.548), whereas those of Jesson and Madden, 
which referred to liquid Al in contact with the (100) face 
of its solid fee phase, gave A w 2.1 A (which is close to 
the interlayer distance in the solid) and therefore leads 
to a cr/A w 0.613. This latter ratio is close to the range 
obtained for the alkalis and therefore it foretells an op- 
posite phase between the longitudinal ionic and valence 
electronic DPs, in accordance with Jesson and Madden's 
results.^ 
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FIG. 17: Valence electronic density profiles obtained by 
superposition of model densities with different a/A values, 
namely (j/A= 0.45 (full line), 0.575 (broken line) and 0.625 
(dotted line) . The thick line represents the longitudinal ionic 
density profile and the inset for K depicts the two outermost 
maxima. 



Experimental study of the LV interface is usually per- 
formed by X-ray reflectivity and/or grazing incidence X- 
ray diffraction techniques. Both methods measure the 
total electronic density profile which is the sum of both 
the core and valence electronic density contributions. 
Whereas the valence electronic DP is obtained by the OF- 
AIMD method, the corresponding core electronic DP has 
been calculated as follows. Using the OF-AIMD gener- 
ated ionic positions we have placed at each ionic position 
the core electronic density (already computed in the pro- 
cess of calculating of the local pseudopotential) and we 
have carried out a configurational average. Once again. 
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FIG. 18: Valence (dotted line), core (dashed line) and total 
(fuU line) electronic density profiles (relative to the valence 
bulk value) for liquid Mg and Si. The longitudinal ionic DP 
is indistinguishable from the total electronic one. 



we note that the core electronic density is a spherically 
symmetric function which has previously been integrated 
over its x and y-dependencies yielding a funcion nc{z), 
which is the one effectively used in the calculation of the 
core electronic DP. Figure ITSl shows, for the Mg and Si 
slabs, a comparison among the obtained valence, core 
and total electronic DPs. Since the core electronic den- 
sities ric^z), are rather narrow funcions, their superposi- 
tion leads to a core electronic DP which stands clearly 
in phase with the longitudinal ionic DP. Consequently, 
when it is added to the valence electronic DP it yields 
a total electronic DP whose phase-shift with respect to 
the longitudinal ionic DP depends on both the relative 
weight of the core electronic DP (always in phase) and the 
valence electronic DP (any phase is possible) as well as 
on the amplitude of the oscillations in both DPs. For the 
particular cases of Mg and Si, which are plotted in Figure 
1181 their valence electronic DPs have relative weights of 
1/6 and 2/7 respectively and, moreover, the oscillations 
in their core electronic DPs are of substantially larger am- 
plitude than those in their respective valence electronic 
DPs. Both features cooperate so as to produce a total 



electronic DP which practically coincides with longitudi- 
nal ionic DP. Indeed, this is a common trait in all the 
metallic systems considered in this work. Even liquid Li, 
which a priori could be more prone to changes because 
its valence electronic DP has a relative weight of 1/3 and 
also stands in opposite phase with the core electronic DP, 
has a total electronic DP which closely matches the ionic 
one. 



IV. CONCLUSIONS 

Results have been reported for the structure of the 
LV interface in several simple sp-bonded liquid metals. 
They have been obtained by an ab initio molecular dy- 
namics method based on the density functional theory. 
Although it employs an approximate electronic kinetic 
energy functional, the large variations in electron density 
associated with the liquid- vapour interface are accounted 
for selfconsistently in the forces acting among the ions. 
The MD simulations were performed by using simulation 
slabs composed of 2000 ions, which are wide enough to 
discard possible interference effects between the two free 
surfaces. 

The calculated longitudinal ionic and electronic DPs 
exhibit clear oscillations, with the ionic ones lasting for 
at least four layers into the bulk liquid. The wavelength 
of the ionic oscillations shows good scaling with the radii 
of the associated Wigner-Seitz spheres; conversely no def- 
inite relationship with electronic parameters has been 
found. Furthermore, the metals with the greater valence 
display ionic DPs with the higher amplitude of the outer 
oscillation as well as a steeper decaying tail in the LV 
interface region. 

Several structural properties were calculated, with spe- 
cial attention devoted towards its evolution along the 
outermost layers. Results were presented for the trans- 
verse pair distribution functions which in turn are used 
to analyze the atomic rearrangements occurring at LV 
interface. The coordination numbers remain practically 
constant for most of the slab and only very close to the 
LV interface there is a reduction related to the absence 
of neighbors outside the interface. However, when com- 
pared properly, we have found that the angular distribu- 
tion of the nearest neighbors doesn't change significantly 
across the interface. Moreover, it is shown that the inter- 
face induces weak structural rearrangements as compared 
with and ideally, step-like, terminated surface. 

The valence electronic DPs show oscillations near the 
LV interface which are much weaker than those of the 
associated ionic DP. Moreover, it was found that the va- 
lence electronic DP is practically reproduced by super- 
posing, at the ionic sites, the pseudoatomic valence den- 
sities calculated in the process of constructing the local 
pseudopotential. This suggests that for each system, the 
main features of its self-consistent valence electronic DP 
are already embodied in the characteristics of the corre- 
sponding pseudoatom valence density. 
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We have also analized the mechanisms behind the rel- 
ative phases of the oscillations in the ionic and valence 
electronic DPs. It is found that those phases evolve 
from opposite phase in the alkalis to almost in-phase for 
Si. This is in stark contrast with the accepted wisdom, 
namely that the electronic and ionic DPs should oscillate 
in opposite phase. An explanation is provided in terms of 
the size of the pseudoatoms (ion plus valence electronic 
cloud) relative to the distance between consecutive layers 
of the ionic profile. 
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